{
 "metadata": {
  "language": "Julia",
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Nonlinear mixed-effects models in Julia"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Linear mixed-effects models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A linear mixed-effects model is characterized by the distribution of\n",
      "two vector-valued random variables: the $n$-dimensional response, $Y$,\n",
      "and the $q$-dimensional random effects vector, $B$.  The unconditional\n",
      "distribution of $B$ is multivariate normal\n",
      "$$\n",
      "    B\\sim N(0,\\Sigma_\\theta)\n",
      "$$\n",
      "as is the conditional distribution of $Y$ given $B=b$\n",
      "$$\n",
      "    (Y|B=b)\\sim N(X\\beta+Zb, \\sigma^2 I)\n",
      "$$\n",
      "\n",
      "In the [MixedModels](https://github.com/dmbates/MixedModels) package for \n",
      "[Julia](http://julialang.org) we represent the covariance matrix in the unconditional distribution of $B$ as\n",
      "$$\n",
      "    \\Sigma_\\theta=\\sigma^2\\Lambda_\\theta\\Lambda_\\theta'\n",
      "$$\n",
      "where $\\Lambda_\\theta$ is the $q\\times q$ _relative covariance factor_.\n",
      "\n",
      "For given values of $\\theta$ and $\\beta$ we solve a penalized linear least squares problem of the form\n",
      "$$\n",
      "    r^2_{\\beta,\\theta}=\\min_u \\|(y-X\\beta)-Z\\Lambda_\\theta u\\|^2 + \\|u\\|^2\n",
      "$$\n",
      "for which we compute the sparse Cholesky factorization\n",
      "$$\n",
      "    L_\\theta L_\\theta' = \\Lambda_\\theta'Z'Z\\Lambda_\\theta + I\n",
      "$$\n",
      "where $L_\\theta$ is a lower-triangular sparse matrix.\n",
      "Because $L_\\theta$ is triangular, we can easily evaluate its determinant as the product of its diagonal elements.\n",
      "By construction these diagonal elements are positive and the log-determinant, $\\log(|\\Lambda_\\theta'Z'Z\\Lambda_\\theta + I|)=2\\log(|L_\\theta|)$ is easily evaluated.\n",
      "\n",
      "The log-likelihood for a linear mixed-effects model, $\\ell(\\beta,\\theta|y)$, can be expressed as\n",
      "$$\n",
      "    -2\\ell(\\beta,\\theta|y) = \\log(|\\Lambda_\\theta'Z'Z\\Lambda_\\theta + I|)+n\\left[1+\\log\\left(\\frac{2\\pi r^2_{\\beta,\\theta}}{n}\\right)\\right]\n",
      "$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Nonlinear mixed-effects models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The formulation of a nonlinear mixed-effects model (NLMM) is nearly the same as that of an LMM except that the mean of the conditional distribution $Y|B=b$ is a nonlinear function of the parameters $\\beta$ and the random effects, $b$ or, equivalently, the _spherical_ random effects $u$ where $b=\\Lambda_\\theta u$.\n",
      "\n",
      "The nonlinear model function, $f$, is a function of a $k$-dimensional model parameter, $\\phi$ and covariates. In our formulation we write the `vec` of the $n\\times k$ _model parameter matrix_, $\\Phi$ as a linear predictor\n",
      "$$\n",
      "    vec(\\Phi)=X\\beta+Z\\Lambda_\\theta u\n",
      "$$\n",
      "and evaluate\n",
      "$$\n",
      "    \\mu_{Y|U=u}=f(T,\\Phi)\n",
      "$$\n",
      "where $T$ is a matrix of covariate values, _time_ in the case of a population pharmacokinetic model.\n",
      "\n",
      "The penalized linear least squares (PLS) problem solved for each evaluation of the objective in the case of a linear mixed-effects model is replaced by a penalized nonlinear least squares (PNLS) problem\n",
      "$$\n",
      "    r^2_{\\beta,\\theta}=\\min_u \\|y-\\mu_{Y|U=u}\\|^2+\\|u\\|^2\n",
      "$$\n",
      "in an NLMM.  \n",
      "The optimizer, $\\tilde{u}=\\arg\\min_u r^2_{\\beta,\\theta}$, is the _mode_ of the distribution of $U$ given $Y=y$ and the current values of the parameters, $\\beta$ and $\\theta$.\n",
      "\n",
      "In general there will not be a closed-form expression for the log-likelihood.  However, the _Laplace approximation_ to the log-likelihood is exactly the same expression as the log-likelihood for an LMM. The Laplace approximation is equivalent to a one-point adaptive Gauss-Hermite approximation to the integral that defines the log-likelihood.  This approximation is _adaptive_ in the sense that the Gauss-Hermite approximation is taken at the mode of the conditional distribution, not at the mean of the unconditional distribution."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A simple example - first-order kinetics, 1 compartment, single bolus unit dose"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One of the simplest pharmacokinetics models is that for a one compartment model with a single bolus dose.\n",
      "The predicted concentration at time $t$ is given by\n",
      "$$\n",
      "        c(t)=V\\exp(-Kt)\n",
      "$$\n",
      "where V is the volume of distribution and K is the elimination rate constant.\n",
      "We represent an instance of such a model with the user-defined type `logsd1`."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Fitting the model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Simulated data from a one compartment model with a single bolus dose are provided in a compressed `csv` file in the `data` directory of the `NLreg` package.  We read these data as a `DataFrame` object.  The `head` function shows us the first few rows."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "using NLreg, DataFrames\n",
      "sd1 = readtable(Pkg.dir(\"NLreg\",\"data\",\"sd1.csv.gz\"));"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We create a nonlinear regression model from these data as"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = BolusSD1(:(CONC ~ TIME), sd1);\n",
      "names(m)'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "1x5 Array{Any,2}:\n",
        " :t  :y  :mu  :resid  :tgrad"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "Ignoring the between-subject differences we can fit the model as"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nl = NonlinearLS(m)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Residual sum of squares at estimates: 110"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "Model fit by nonlinear least squares to 580 observations\n",
        "2x4 DataFrame:\n",
        "        parameter estimate   stderr t_value\n",
        "[1,]          \"V\"  1.14295 0.049565 23.0596\n",
        "[2,]          \"K\" 0.245682 0.020241 12.1378\n",
        "\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".597\n",
        "Residual standard error = 0.43743 on 578 degrees of freedom\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `SimpleNLMM` type represents a nonlinear mixed-effects model with a single grouping factor for the random effects (`ID` in this case), and random-effects for each of the nonlinear model parameters. The third and fourth arguments to the external constructor are the initial values for the template $\\Lambda$ matrix and for the fixed-effects parameters $\\beta$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp = SimpleNLMM(nl,vector(sd1[:ID]),Diagonal(ones(2)));\n",
      "prss!(pp,0.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "110.59728633992346"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `prss!` function evaluates the penalized residual sum of squares, $r^2_{\\beta,\\theta}(u)$ at the initial parameter estimates.  The fixed-effects parameter, $\\beta$, is the parameter vector from the nonlinear least squares fit, the template for the relative covariance factor, $\\Lambda$, is initialized to the identity and the random-effects, $u$, are initialized to zero."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit(pp;verbose=true);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "LoadError",
       "evalue": "fit not defined\nat In[4]:1",
       "output_type": "pyerr",
       "traceback": [
        "fit not defined\nat In[4]:1"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The individual values of the nonlinear model parameters, $\\phi_i,i=1,\\dots,200$ are stored as the `phi` element of the model structure."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp.phi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "2x200 Array{Float64,2}:\n",
        " 0.549264  0.907821  0.897507  0.678659  \u2026  1.27107   0.489555  0.914113\n",
        " 0.263226  0.267552  0.316919  0.330285     0.191895  0.325349  0.331993"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "names(pp)'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "1x14 Array{Any,2}:\n",
        " :m  :inds  :nrep  :lambda  :L  :beta  \u2026  :nAGQ  :mxpnls  :minfac  :tolsqr"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp.lambda"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "2x2 Diagonal{Float64}:\n",
        "diag:  4.55744  0.748434"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pp.beta'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "1x2 Array{Float64,2}:\n",
        " 1.19333  0.294698"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}